#-------------------------------------------------------------------
#this script include the definition of partial credit model/generalized partial credit model
#-------------------------------------------------------------------

# -----------------------------------------------------------------------
#   P_gpcm computes the probability of each response for an item
#   that is well fit by the generalized partial credit model.  The program
#   has as input a vector of ability parameters, theta, the discrimination
#   parameter for the item, a, the dificulty parameter for the item, b, and
#   the threshold parameters in the vector d.  The output is a matrix of  
#   probabilities for each response for each ability in theta.
#------------------------------------------------------------------------

P_gpcm <- function(theta,a,b,d){

	#   Determine number of abilities.
	n=length(theta);

	#   Determine the number of score categories.
	m=length(d);

	#   Determine values of exponents
	num = matrix(0.,ncol=m,nrow=n)
	for(i in 1:m){
    		s=0;
    		for(k in 1:i){
        		s=s+1.7*a*(theta-b+d[k]);
    		}
    	num[,i]=exp(s);
	}

	denom=as.matrix(rowSums(num));
	ones=matrix(1,nrow=1,ncol=m)
	den=denom%*%ones
	p=num/den;
}


#--------------------------------------------------------------------------------
# iteminfo_gpcm calculate the item info for generalized partial credit model.
#--------------------------------------------------------------------------------
iteminfo_gpcm<-function(theta,a,b,d){  
        n=length(theta);
        m=length(d);
	Nitem = length(a);
        D = 1.7;
        kk = seq(1,m);
        one = matrix(1,nrow=n,ncol=1);
        P=P_gpcm(theta,a,b,d);
	  Eut = P%*%kk # expectation
        s=0;
        for(k in 1:m){
           s = s + (k*one-Eut)^2*P[,k];
	  }
        inf=D^2*a^2*s;
}


#--------------------------------------------------------------------------------
# testinfo_gpcm calculate the test info for generalized partial credit model.
#--------------------------------------------------------------------------------
testinfo_gpcm<-function(theta,a,b,d){
      Nitem = length(a);
      s=0;
      for(i in 1:Nitem){
         s = s + iteminfo_gpcm(theta,a[i],b[i],d[i,]);
      }
      inf = s ;
}

#-----------------------------------------------------------------
#standard error of gpcm
#-----------------------------------------------------------------
sderr_gpcm<-function(theta,a,b,d){
      res=1./sqrt(testinfo_gpcm(theta,a,b,d));
}



# -----------------------------------------------------------------------
#   P_2ppcm computes the probability of each response for an item
#   that is well fit by the two parameter partial credit model.  The program
#   has as input a vector of ability parameters, theta, the discrimination
#   parameter for the item, a, and
#   the threshold parameters in the vector d.  The output is a matrix of  
#   probabilities for each response for each ability in theta.
#   the input 
#------------------------------------------------------------------------

P_2ppcm <- function(theta,a,d){

	#   Determine number of abilities.
	n=length(theta);

	#   Determine the number of score categories.
	m=length(d);

	#   Determine values of exponents
	num = matrix(0.,ncol=m,nrow=n)
	for(i in 1:m){
    		s=0;
    		for(k in 1:i){
        		s=s+1.7*a*(theta+d[k]);
    		}
    	num[,i]=exp(s);
	}

	denom=as.matrix(rowSums(num));
	ones=matrix(1,nrow=1,ncol=m)
	den=denom%*%ones
	p=num/den;
}


#-------------------------------------------------------------------------------
# iteminfo_2ppcm calculate the item info for 2 parameter credit model.
# input: theta as a vector/scalar; a as scalar, d as vector, see testinfo_2ppcm
#--------------------------------------------------------------------------------
iteminfo_2ppcm<-function(theta,a,d){  
        n=length(theta);
        m=length(d);
	Nitem = length(a);
        D = 1.7;
        kk = seq(1,m);
        one = matrix(1,nrow=n,ncol=1);
        P=P_2ppcm(theta,a,d);
	Eut = P%*%kk # expectation
        s=0;
        for(k in 1:m){
           s = s + (k*one-Eut)^2*P[,k];
	  }
        inf=D^2*a^2*s;
}


#--------------------------------------------------------------------------------
# testinfo_2ppcm calculate the test info for 2 parameter partial credit model.
#--------------------------------------------------------------------------------
testinfo_2ppcm<-function(theta,a,d){
      Nitem = length(a);
      s=0;
      for(i in 1:Nitem){
         s = s + iteminfo_2ppcm(theta,a[i],d[i,]);
      }
      inf = s ;
}



#---------------------------------------------------------
#calculate the mean P of 2 parameter partial credit model-
#input is the theta, frequency of each theta, a and d    - 
#---------------------------------------------------------

P_mean_OE<-function(theta,freq,a,d){
      res=matrix(0.,nrow=dim(d)[1],ncol=dim(d)[2])
	for(j in 1:length(a)){
               s=0.;
               for(i in 1:length(theta)){
                        s=s+P_2ppcm(theta[i],a[j],d[j,])*freq[i]
		} 
	res[j,] = s/sum(freq)
	}
	return(res)
}



#-------------------------------------------------------
# calculate the TCC for 2 parameter partial credit model
#-------------------------------------------------------

tcc_2ppcm<-function(theta,a,d){
   Ntheta=length(theta)
   Nitem=length(a)
   tcc=matrix(0,nrow=Ntheta,ncol=dim(d)[2])
   for(i in 1:Ntheta){
     s=0.;
     for(j in 1:Nitem){
        s=s+P_2ppcm(theta[i],a[j],d[j,])
      }
      tcc[i,]=s   
   }
   
   return(tcc);
}


#-----------------------------------------------------------------
#standard error of 2ppcm
#-----------------------------------------------------------------
sderr_2ppcm<-function(theta,a,b,d){
      res=1./sqrt(testinfo_2ppcm(theta,a,b,d));
}
